Algorithms for multivariate image analysis and other large-scale applications of multivariate curve resolution (MCR) typically employ iterative, linearly constrained alternating least squares (ALS) algorithms in their solution. The efficiency with which the original problem can be solved, therefore, rests heavily on the underlying constrained least squares solver. Efficiency of solution is particularly important with large-scale MCR applications. In this context, the term “large-scale” can be used to denote problems in which the number of observation vectors is much greater than the number of variables to be estimated from each observation. A typical spectral imaging application, for instance, might require that a complete spectrum be measured at each pixel in a spatial array of a large number (e.g., of order 106) of total pixels in order to estimate the concentrations of a modest number (e.g., of order 10) of chemical components at each pixel. A least squares problem having multiple observation vectors is denoted herein as a multiple right hand side (RHS) problem. See, e.g., P. Kotula and M. Keenan, “Automated analysis of SEM X-ray spectral images: A powerful new microanalysis tool,” Microsc. Microanal. 9, 1 (2003); U.S. Pat. No. 6,584,413 to Keenan and Kotula; and U.S. Pat. No. 6,675,106 to Keenan and Kotula. Descriptions of the current state-of-the-art for solving constrained least squares problems in chemical applications can be found in the literature. See, e.g., C. L. Lawson and R. J. Hanson, Solving Least Square Problems, Prentice-Hall, Inc., Englewood Cliffs, N.J. (1974); A. Bjorck, Numerical Methods for Least Squares Problems, SIAM, Philadelphia (1995); R. Bro and S. DeJong, “A fast non-negativity-constrained least squares algorithm,” J. Chemometrics 11, 393 (1997); and M. Van Benthem et al., “Application of equality constraints on variables during alternating least squares procedures,” J. Chemometrics 16, 613 (2002).
Given its fundamental role in solving a variety of linearly constrained optimization problems, much attention has been paid to methods for estimating least squares models subject to non-negativity constraints in particular. The method of non-negativity constrained least squares (NNLS) is an important tool in the chemometrician's toolbox. Often, it is used in ALS algorithms when the variables to be estimated are known to be non-negative (such as fluorescence spectra or chemical concentrations). The simplest way to implement these non-negativity constraints is to solve the corresponding unconstrained least squares problem and then overwrite any negative values with zeros. Unfortunately, this overwriting method does not employ an appropriate least squares criterion so one is left with a solution that has ill-defined mathematical properties. ALS algorithms based on the overwriting method, for instance, do not necessarily lower the cost function on successive iterations and, therefore, are not guaranteed to converge to a least squares minimum. Fortunately, there exist methods that solve the NNLS problems in a mathematically rigorous fashion. These methods impose non-negativity on the solution while minimizing the sum of squared residuals between the data being fit and their estimates in a true least-squares sense.
While rigorous methods ensure true least squares solutions that satisfy the non-negativity constraints, they can be computationally slow and demanding. In addition, rigorous methods are more difficult to program than the simple overwriting method. These are probably the principal reasons many practitioners of techniques that employ NNLS choose not to use rigorous methods. The standard NNLS algorithm originally presented by Lawson and Hanson was designed to solve the NNLS problem for the case of a single vector of observations. See Lawson and Hanson, Chapter 23. When applied in a straightforward manner to large-scale, multiple RHS problems, however, performance of the standard NNLS algorithm is found to be unacceptably slow owing to the need to perform a full pseudoinverse calculation for each observation vector. Bro made a substantial speed improvement to the standard NNLS algorithm for the multiple RHS case. Bro's fundamental realization was that large parts of the pseudoinverse could be computed once but used repeatedly. Specifically, his algorithm precomputes the cross-product matrices that appear in the normal equations formulation of the least squares problem. Bro also observed that, during ALS procedures, solutions tend to change only slightly from iteration to iteration. In an extension to his NNLS algorithm that he characterized as being for “advanced users,” Bro retained information about the previous iteration's solution and was able to extract further performance improvements in ALS applications that employ NNLS. These innovations led to a substantial performance increase when analyzing large multivariate, multi-way data sets. However, even with Bro's improvements, a need remains for further performance improvements and for a fast algorithm that has general applicability to the solution of linearly constrained, multiple RHS least squares problems.
The present invention is directed to a fast combinatorial algorithm that has improved performance compared to prior algorithms for solving linearly constrained least squares problems. The combinatorial algorithm exploits the structure that exists in large-scale problems in order to minimize the number of arithmetic operations required to obtain a solution. Based on combinatorial reasoning, this algorithm rearranges the calculations typical of standard solution methods. This rearrangement serves to reduce substantially the amount of computation required in a typical MCR problem. The combinatorial algorithm has general applicability for solving least squares problems subject to a wide variety of equality and inequality constraints.
The application of the combinatorial algorithm to the solution of the NNLS problem is described herein in detail. The combinatorial algorithm has been tested by performing calculations that are typical of ALS applications on spectral image data. The spectral image data sets varied widely in size and, in all cases, the combinatorial algorithm outperforms both Bro's simple and advanced-user modes, and it does so without the attendant overhead required in the latter case. Using the overwriting method as a point of reference, the combinatorial algorithm also exhibits much better scaling behavior for problems with very large (e.g., greater than 105) numbers of RHS vectors.